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ABSTRACT 

Stellar limb-darkening impacts a wide range of astronomical measurements. The accu¬ 
racy to which it is modelled limits the accuracy in any covariant parameters of interest, 
such as the radius of a transiting planet. With the ever growing availability of precise 
observations and the importance of robust estimates of astrophysical parameters, an 
emerging trend has been to freely fit the limb-darkening coefficients (LDCs) describing 
a limb-darkening law of choice, in order to propagate our ignorance of the true intensity 
profile. In practice, this approach has been limited to two-parameter limb-darkening 
laws, such as the quadratic law, due to the relative ease of sampling the physically 
allowed range of LDCs. Here, we provide a highly efficient method for sampling LDCs 
describing a more accurate three-parameter non-linear law. We first derive analytic 
criteria which can quickly test if a set of LDCs are physical, although naive sampling 
with these criteria leads to an acceptance rate less than 1%. We then show that the 
loci of allowed LDCs can be transformed into a cone-like volume, from which we are 
able to draw uniform samples. We show that samples drawn uniformly from the conal 
region are physically valid in 97.3% of realizations and encompass 94.4% of the volume 
of allowed parameter space. We provide python and fortran code (ldc3) to sample 
from this region (and perform the reverse calculation) at this http, which also includes 
a subroutine to efficiently test whether a sample is physically valid or not. 
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1 INTRODUCTION 

Stellar limb-darkening affects a wide range of astronomical 
observations, such as exoplanetary transits (Mandel & Agol 
2002), microlensed light curves (e.g. Witt 1995; Zub et al. 
2011), rotational modulations (Kipping 2012), ellipsoidal 
variations (Morris & Naftilan 1993), interferometric images 
of stars (e.g. Aufdenberg et al. 2005) and eclipsing bina¬ 
ries (Kopal 1950). When interpreting such observations, 
the assumed shape of the limb-darkening intensity profile 
may significantly affect the inferred parameters of interest 
(Csizmadia et al. 2013). Consequently, an accurate treat¬ 
ment of limb-darkening is crucial, even when the effect itself 
is a “nuisance” phenomenon. 

In many practical cases, limb-darkening is treated by de¬ 
scribing the stellar intensity profile with a closed-form ana¬ 
lytic model. This model is usually called the “limb-darkening 
law”, which is designed to provide an accurate analytic ap¬ 
proximation the true profile. The major advantage of mod¬ 
elling the intensity profile analytically is that the astronom¬ 
ical phenomena under investigation may also be modelled 
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analytically, offering computational expedience and greater 
physical insight. As an example, in the field of exoplanet 
transits, the shape of a transit light curve can be described 
with a closed-form, analytic model under the assumption of 
a polynomial-based description of the stellar intensity profile 
(Mandel & Agol 2002; Gimenez 2006). 

All of the commonly used limb-darkening laws may be 
described as a linear sum of one or more simple functions, 
each of which is weighted by a so-called limb-darkening co¬ 
efficient (LDC). For example, the popular quadratic limb- 
darkening law describes the intensity profile as a quadratic 
series expansion with respect to the angle between the 
line of sight and the emergent intensity (fl) and has two 
LDCs (Kopal 1950). These LDCs control the shape of the 
intensity profile, subject to the flexibility granted by the 
limb-darkening law’s complexity. It is therefore necessary 
to make at least two major decisions with how to treat 
limb-darkening; what limb-darkening law should be used and 
what LDCs should be assigned to this law? 

A typical approach for assigning LDCs is to adopt a 
set of so-called “theoretical” LDCs. In this framework, one 
first simulates an intensity profile using a sophisticated stel¬ 
lar atmosphere model at a particular wavelength or over a 
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chosen integrated bandpass. One then takes this simulation 
and regresses the limb-darkening law of choice to it by find¬ 
ing the maximum likelihood set of LDCs (in some instances 
the LDCs are restricted to ensure physically sound profiles). 
Since the simulated profile is sensitive to several parame¬ 
ters defining the stellar surface (e.g. effective temperature, 
metallicity, surface gravity), several groups have produced 
grid tabulations of LDCs for a wide range of plausible in¬ 
puts and assumed limb-darkening law (e.g. see Van Hamme 
1993, Dfaz-Cordovez et al. 1995, Claret 2000 and Sing 2010). 

In recent years, there has been a shift away from using 
theoretical LDCs in favour of freely htting LDCs simulta¬ 
neous to the exploration of the other parameters of inter¬ 
est (e.g. Knutson et al. 2007, Kipping & Bakos 201 la,b and 
Kreidberg et al. 2014). This approach, enabled by advances 
in modern computers, allows one to propagate our ignorance 
regarding the shape of the intensity proHle into the estima¬ 
tion of the parameters of interest, leading to more accurate 
estimates parameter uncertainties. Such a procedure also al¬ 
lows one to decouple from any possible errors in the stel¬ 
lar atmosphere models themselves, which have been shown 
to sometimes be inconsistent with observed limb-darkening 
profiles (Howarth 2011; Epinoza & Jordan 2015). 

One major challenge when attempting to freely ht LDCs 
is that many combinations of the LDCs are unphysical. For 
example, a specific choice of LDCs may result in the inten¬ 
sity profile being occasionally negative. These regions of un¬ 
physical parameter space may coincide with likelihood min¬ 
ima, leading to erroneous parameter posterior distributions. 
Softer, extended likelihood minima (typical of low signal- 
to-noise data) may spill over into the unphysical parameter 
space, causing the final parameter posteriors to be marginal¬ 
ized over large swaths of unphysical parameter space leading 
to unnecessarily swollen error bars. 

There are two solutions to this problem. The first is 
to define a set of criteria which quickly allow us to iden¬ 
tify forbidden combinations of the LDCs. Armed with such 
criteria, one could then test each realization and accept or 
reject the realization accordingly. This works fine unless the 
volume of forbidden parameter space is large, in which case 
one will severely impede a regression algorithm’s efficiency 
since most trials are being rejected as unphysical. A more 
elegant solution is to directly sample exclusively (or at least 
efficiently) from the physically allowed volume of parame¬ 
ter space (and also in such a way that the LDCs can be 
uniformly, or nearly uniformly, sampled). This provides the 
benehts of a very high efficiency, physically sound priors and 
no need for testing criteria at each realization. However, this 
approach comes at the cost of the initial intellectual invest¬ 
ment to actually solve how to perform efficient sampling in 
the hrst place. 

Because direct sampling of physical LDCs is challeng¬ 
ing, the most complex limb-darkening law for which this feat 
has yet been achieved is the simple quadratic law, where the 
intensity profile of the star is given by 

= ( 1 ) 

where /(I) is the specihc intensity at the centre of the 
disc. III and U 2 are the quadratic LDCs and g is the cosine of 
the angle between the line of sight and the emergent inten¬ 
sity. In Kipping (2013), we showed that physically sound ui 


and U 2 LDCs reside within a triangle on the ui-U 2 plane, 
from which uniform samples can be easily drawn by re¬ 
parameterizing to qi = (ill + 112 )^ and q2 = 0.5mi(mi + 112 )^^. 
Other simple two-parameter laws were considered as well in 
this work. 

Despite these successes with two-parameter laws, these 
laws are fundamentally limited in their ability to accurately 
describe realistic limb-darkening prohles. This translates to 
systematic uncertainty in any and all model parameters 
which are covariant with the limb-darkening properties. This 
may, for example, limit our ability to accurately measure the 
radius of a transiting exoplanet and thus make inferences 
of its composition. We are therefore motivated to extend 
the Kipping (2013) analysis to a more sophisticated limb- 
darkening law. 

The most accurate closed-form limb-darkening law is 
the Claret (2000) four-parameter non-linear law, described 
by 

/(g)//(l) = l-ci(l-g^/^)-c 2 (l-g) 

-C3(l-/t3/2)_c4(l_^2). (2) 

In numerous independent studies, this law has been 
found to provide the most accurate description of simu¬ 
lated intensity profiles (e.g. Claret 2000, Sing 2010 and 
Magic et al. 2015) versus competing models. This is, how¬ 
ever, not surprising since this law also utilizes the great¬ 
est number of LDCs. Analytically identifying the unphysi¬ 
cal combinations of these four LDCs remains an outstanding 
and formidable challenge. A more tractable problem that we 
consider in this work is the allowed volume of LDCs (and 
methods to directly sample from said volume) in the case of 
the Sing et al. (2009) three-parameter law: 

/(g)//(l) = l-c2(l-g)-c3(l-g^/2)-c4(l-g^). (3) 

Sing (2010) argues that dropping the cj term is mo¬ 
tivated by solar data (Neckel & Labs 1994) and 3D stel¬ 
lar models (Bigot et al. 2006), which show that /(g) varies 
smoothly at small g, meaning that a g^^^ term is superflu¬ 
ous. We therefore argue that the Sing et al. (2009) law offers 
both a significant improvement in accuracy and yet is simple 
enough for us to analytically constrain the allowed LDCs. 

2 PHYSICAL CRITERIA 
2.1 Physical Conditions 

We define the following physical conditions with respect to 
a limb darkened stellar intensity proHle: 

(I) an everywhere-positive intensity profile, 

(II) a monotonically decreasing intensity profile from the 
centre of the star to the limb, 

(III) the intensity prohle has a negative curl at the limb. 

Physical conditions I and II are the same two imposed 
in our previous paper. Kipping (2013). As noted in that 
work, limb brightening is possible for narrow-band observa¬ 
tions (e.g. see Schlawin et al. 2010) and such behaviour is 
not considered in this work either. Physical condition III 
is motivated by the expectation that the intensity rapidly 
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drops off towards the limb (and is discussed in more de¬ 
tail later in §2.5). Throughout this work, we refer to a set of 
LDCs satisfying these three conditions as being physical, and 
LDCs otherwise are defined as unphysical. In what follows, 
we explore the consequences of these simple constraints. 


lim/(/i) < lim/(/i). (11) 

/i->0 

This condition, derived by physical condition II, pro¬ 
vides our fourth criterion. 


2.2 Criterion A 

Physical condition I demands that I{g) >0V0<^<1. We 
begin by evaluating this condition at two extrema of /t —S’ 1 
and /t —> 0, in a similar manner to the approach adopted in 
Kipping (2013): 

lim / = 1 > 0, 

lim / = 1 — C 2 — C 3 — C 4 > 0. (4) 

The upper line clearly has no constraining power, but 
the second line provides our first criterion of 

C2 C'i C4 1. (5) 

2.3 Criteria B and C 


C2 -f C3 -t- C4 > 0. (12) 

2.5 Criterion E 

Consider the behaviour of the intensity profile at the limb. 
By virtue of condition II, the intensity profile must be de¬ 
creasing as we approach the boundary. We therefore expect 
a negative gradient with respect to r, or equivalently a pos¬ 
itive gradient with respect to fl, since dr/dfi is always neg¬ 
ative. We also expect that at the limb the gradient of the 
gradient (i.e. the curl) is negative. This is consistent with 
the asymptotic-like behaviour expected due to foreshorten¬ 
ing near the limb, causing the gradient to become ever-more 
negative and defines physical condition III. Note that a neg¬ 
ative curl with respect to r is equivalent to a negative curl 
with respect to fl, since we now multiply by (dr/dfl) twice, 
leading to a double negative. The curl may be expressed as 


Next, we consider physical condition II, which demands that 
dl/dfl > 0 y 0 < fl < i: 

= C 2 + ^C 2 fl^/^+ 2 c 4 fl. (6) 

As was done in the previous subsection, let us evaluate 
the above in the extreme cases of /i —>■ 1 and fl^O, yielding 


dfl^ 


^C3+2c4fl^E^ 


At the limb then (/i —> 0), we expect that 


lim 



< 0 , 


which defines criterion E, 


(13) 


(14) 


,■ dl{fi) 3 „ 

hm ^— = C2 + -C3 -|-2 c4 > 0 , 
g->l dfl 2 


lim 

g^O 


dljp) 

dfl 


= C2 > 0. 


C3 < 0. 


(7) 


2.6 Criterion F 


(15) 


These two expressions provide our criteria B and C, 
which are respectively given by 

2c2 T 3^3 T 4c4 >0, (8) 

and 


Physical condition I requires that I{fl) is everywhere pos¬ 
itive. Combining I with II implies that I(fl) must be less 
than unity everywhere, which one may consider to be phys¬ 
ical condition I’. Writing this out along with II, one may 
show that 


C2 > 0. (9) 

2.4 Criterion D 

Physical condition II tells us that the intensity decreases 
from the centre of the star to the limb. This implies that the 
intensity everywhere (except at the centre of the star) is less 
than that present at the centre of the star, or mathematically 
that 

I{fl)<\iml{fi). (10) 

One simple closed-form result from this constraint oc¬ 
curs by comparing the intensity at the limb to the centre 
via: 


-C2fl-C3,fl'^^^-C4fl^ > -C2-C3-C4, (16) 

2 4 

-C2fl+C3fl^i^ + -C4fl^ >0. (17) 

Adding the two inequalities shown above cancels out 
the c^fl^E terms and leaves us with a quadratic equation: 

C4fl^-C2fl>-3 {c2 + C2 + C4). (18) 

From Criterion A, we know that the sum of the coeffi¬ 
cients must be less than unity, implying that in the limit of 
fl ^ 1, we have, 

C4 — C2 > —3. (19) 

Starting from criterion B and invoking criterion E, we 
can also show that: 
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2c 2 “t" 3 c3 -\- 4c4 > 0, 

3c3 > — 2 c 2 — 4c4, 

— 2 c 2 — 4c4 < 3c3 < 0 , 

C 2 + 2 c 4 > 0 . ( 20 ) 

Summing Equations 19 & 20 together yields 


Cl > 


32c4 


(27) 


2.8 Summary of Analytic Criteria 

To summarize, our seven analytic criteria on the three LDCs 
are 


C 4 >- 1 . ( 21 ) 

Through numerical experimentation, we find that ap¬ 
plying the slightly more conservative bound of C 4 > 0 yields 
a more symmetric loci of allowed points (as discussed later in 
§5), from which it is easier to directly sample. We therefore 
modify criterion F to. 


C 4 > 0. 


( 22 ) 


2.7 Criterion G 

For our final criterion, we begin by considering physical con¬ 
dition II: 


dm 

dp 


> 0 , 


2c4/i -f —-f C2 > 0. 


(23) 


The gradient expressed must be everywhere-positive 
and so let us compute the minimum gradient possible, which 
occurs when the curl equals zero, or when 


(9(2c 4/X +|c3/t'/^ + C2) 

"dp 

^C3/t-'/2 + 2c4=0. (24) 

Therefore, the minimum gradient occurs when p = Pujim 
where we define 


1/2 

m = 


3cj 

8 c 4 


(25) 


In the case where criterion E and F are in effect, then 
C 3 is negative and C 4 is positive meaning that p^i„ is a real 
number. 

The point p„j[„ may or may not be within the range 
0 < p < 1. If indeed it is, then implicitly — 1 < 3 c 3 /( 8 c 4 ) < 0 
and we require that the gradient at this point is positive, 
giving 


if 


3c3 

1 < —4- < 0 then 
8 c 4 


Cl 


32c4 


(26) 


As with criterion F, we find through numerical tests 
in §5 that the loci of points can be made symmetric if we 
impose criterion G under all circumstances, not just when 
— 1 < 3 c3 /(8 c4 ) < 0. We therefore modify criterion G to 


C2 -f C3 4 - C4 < 

2c2 + 3c3 -I- 4c4 > 

Cl > 

C2 -f C3 4 - C4 > 

C 3 < 

C 4 > 

Cl > 

Using the three physical conditions only, we note that 
criterion F should strictly be C 4 > — 1. We have modified 
this criterion to be slightly more conservative so that the 
loci of allowed LDCs can be transformed into a symmetric 
cone shape depicted later in §5. Similarly, criterion G strictly 
only applies when — 1 < 3 c 3 /( 8 c 4 ) < 0 if one uses the three 
physical conditions. We again modify the criterion such that 
it applies under all circumstances, in order to yield a more 
symmetric loci, as shown later in §5. We provide python and 
FORTRAN code (ldc3) to test whether these criteria hold, for 
which the user can also use the unmodified versions of the 
criteria if desired (available at this http). 

In §4, we explore the consequences of these modifica¬ 
tions and perform numerical tests demonstrating the effec¬ 
tiveness of the seven criteria. First though, we calculate the 
allowed maxima/minima on each LDG using the seven cri¬ 
teria, as shown in the next section, §3. 


1, [A] 
0, [B] 
0, [C] 
0, [D] 
0, [E] 
0, [F] 


9c^ 

32c4 


. [G] 


(28) 


3 BOUNDS ON THE LDCs 

In the cases of criteria C, E and F, we have derived a hard 
boundary on the isolated LDCs; for example, C 2 > 0 from 
criterion C. However, it is also useful to know about the 
other extrema; for example, what is the maximum bound 
on C 2 ? Working in the {c 2 ,C 3 ,C 4 } parameter space, which we 
refer to as the “original” reference frame, these bounds define 
the smallest cuboid within which the allowed loci reside. 


3.1 C 3 -C 2 plane 

Starting from the seven criteria, how can we calculate the 
limits on each LDC? We choose to proceed by means of 
numerically-guided analytic reasoning via inspection of the 
projected 2D planes. 

We begin by generating a large number of random 
points drawn uniformly in {C 2 ,C 3 ,C 4 } and then reject any 
point which does not satisfy criteria A-G. The surviving 
points are then saved. We plot the loci of these points on 
the C 3 -C 2 plane in Fig. 1. 

Inspection of Fig. 1 reveals that the loci of points 
are nearly perfectly enveloped by the functions C 3 = 
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Figure 1. Loci of points satisfying criteria A-G, plotted in the 
C 3 -C 2 plane. The pink/blue lines are given by 



C 2 


Figure 2. Loci of points satisfying criteria A-G, plotted in the 
C4-C2 plane. The pink/blue lines are given by 


(4/9)(— 4 c 2 ± y' 2 c 2\/9 — C2), with the exception of a small re¬ 
gion prohibited by criterion E (C 3 < 0). The envelope func¬ 
tion may be derived starting from criteria A and G, as fol¬ 
lows: 

C4 < 1 - C3 - C2, [A] 

9c2 


C 4 > 


32c2 


, [G] 


9ci 

< I-C 3 -C 2 . 

32c2 


(29) 


I to 


which may be re-arranged 

^(^ 3 - ^(-4c2“ V^V9“C2)) 

^ (*^3 - ^-4c2 + a/^a/9-C2)) < 0 

For the above to hold, then C 3 must be larger than one 
of the radicals but always less than the other, giving two 
possible sets of envelope functions. We define the first one 
as: 


(30) 


i/2c2i/9 — C2 = 0, giving C2 = 9. The minimum C3 value can 
be found by minimizing the lower envelope, which occurs at 
C 2 = (3/2)(3-|-2\/2), corresponding to the minimum value of 
C3 = —8 — 63 / 2 , or —16.4853.... This therefore provides two of 
the missing three LDC bounds we seek. 

3.2 C 4 -C 2 plane 

We may repeat this exercise in the C4-C2 plane, and the re¬ 
sulting loci are shown in Fig. 2. As before, two lines appear 
to nearly perfectly envelope the loci of points, which can be 
derived starting from criteria A and G: 

C 3 > (C 2 + C 4 - 1)^, [A] 

„2 ^ 32 C 4 C 2 
^^3 < - Q -> [G] 


(C2-I-C4- 1 )^ < 


9 

32 C 4 C 2 


(33) 


where the last line may now be re-expressed as 


C3 > c 


SjCnv'' 




(^C4- ^(9 + 7c2-4v^^3/9-C2)j 

1 

4c2- 3 /^ 379 - 02 ), 

(31) 

x(c 4 - ^( 9 -b 7 c 2 -|- 43 /^ 3 / 9 -C 2 )j 

1 < 0. (34) 

4c2-|-a/^3/9-C2). 

(32) 

The above requires that C 4 is 

greater than one of the 


In Fig. 1, the pink line depicts c^env(‘' 2 )i whilst the blue 
line depicts c/gjj^(c 2 ). This demonstrates that indeed our 
first guess for the form of the solution is correct. In principal 
though, we note that an alternative solution is c^e„v(c 2 ) < 
C3 < C3/env(c2)- 

From Fig. 1, one can see that there is a minimum al¬ 
lowed C 3 value and a maximum 02 - The maximum C 2 coin¬ 
cides where the two envelope functions meet, meaning that 


radicals but always less than the other. Therefore, we have 
two valid envelope functions and we define the first one as: 


1 


C4 >C4 ,,„^(C2) = -(9-|-7c2-4i/^V 9-C2), 
<^4 < cjenv (C2) = ^ (9 + 7 c2 + 4^^^9^) . 


(35) 

(36) 


These two functions envelope the loci of points simu¬ 
lated earlier. This is apparent from Fig. 2, where the pink 
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-10 -5 

C3 


Figure 3. Loci of points satisfying criteria A-G, plotted in the 
C4-C3 plane. The pink/blue lines are given by 


line denotes c^g„y(c 2 ) and the blue line C 4 env(^ 2 )- We plot the 
^ 4 env('' 2 ) function from the upper-right intersection down to 
hitting the C 4 boundary condition of criterion F. Extending 
the line further back, shown as a dotted line, does not pro¬ 
vide a physical bound on the points, although we note that 
only a very small fraction of points are excluded by doing 
so. In any case, the objective here is merely to determine 
the upper limit on C4. 

We may use the c^j.ny(c 2 ) function to find the maximum 
bound on C4. Maximizing the c^env(‘'2) function by differen¬ 
tiation, we find the curve is maximized at C 2 = 8 , correspond¬ 
ing to C 4 = 9. This provides the last limit needed to define 
the cuboid containing all loci as tightly as possible: 

0 < C 2 < 9, 

—8 —6\/2 < C3 < 0, 

0 < C 4 < 9. (37) 


3.3 C4-C3 plane 

Although we have now derived all of the LDC bounds, for 
the sake of completeness we here consider envelope functions 
bounding the C 4 -C 2 plane, as shown in Fig. 3. 

The pink and blue lines shown in Fig. 3 nearly perfectly 
envelope the loci of allowed points, with the exception of 
criterion E (C 3 < 0) truncating a small fraction of points. As 
before, these functions may be derived from the following: 


C2 < I-C3-C4, [A] 

9ci 


32c4 


< 1 — C3 — C4, 


(38) 


which may be re-expressed as 

32 (c4 - i (4 - 4 c3 - v^^8-16c3-c 2) ) 
x(c4-i(4-4c3-l-v^^8-16c3-c2)) < 0. (39) 

Following the lines of argument used before, this pro¬ 
vides the envelope functions: 

M >C4^^env(c3) = ^(4-4c3 - ^2^8- 16c3-c^), (40) 

M (C3) = ^ (4 - 4c 3 -b 72 ^8 - 16c 3 - c^). (41) 

We note that the two functions meet when 
^ 8 — I 6 C 3 — Cj = 0, occurring at C 3 = —8 4 - 63/2 = 0.4853.... 
However, criterion G truncates these functions, albeit only 
a small fraction of the envelope. 


4 NUMERICAL TESTING OF THE CRITERIA 
4.1 Overview 

Starting from three physically imposed conditions for the 
intensity profile of a limb darkened star, we have derived 
seven criteria which bound the three LDCs parameterizing 
the Sing et al. (2009) limb-darkening law. In this section, we 
test the validity of these criteria in terms of (i) completeness 
and (ii) validity. We define these terms as follows: 

■ Completeness: A fully complete set of criteria require 
that the loci of points for which the physical conditions are 
met also satisfy our seven criteria. 

■ Validity: A fully valid set of criteria requires that that 
the loci of points which meet our seven criteria never break 
the two physical conditions. 

These two tests can be thought of in the following 
way. Incomplete cases imply that our criteria are overly- 
conservative, cropping some parts of physically plausible 
parameter space. Invalid cases imply that our criteria are 
overly-optimistic, erroneously predicting that some parts of 
parameter space are physically plausible. 


4.2 Completeness Tests 

We perform our tests via numerical Monte Carlo simula¬ 
tion. We begin by drawing a uniform random point in the 
parameter space {c 2 ,C 3 ,C 4 }. The expected upper and lower 
bounds on these terms were calculated earlier in §3.2. We 
begin by using these bounds except that we use the orig¬ 
inal C 4 > — 1 constraint rather than the modified criterion 
F. We then take this cuboid and double the lengths of each 
side such that we consider the ranges: —9/2 < C 2 < 27/2, 
- 3 ( 44 - 33 / 2 ) < C3 < ( 4 - 433 / 2 ) and —6 < C4 < 14. These ad¬ 
justments are made to ensure that we explore the full range 
of physically allowed LDCs during this test. 

A sample point is drawn from this cuboid and then 
tested as to whether the physical conditions I, II and III 
are met. In practice, we accomplish this by computing 10^ 
points along the functions I{jl) and /dfl varying g from 
0 to 1 in equal, linear steps. For III we simply test if C 3 < 0, 
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since this condition only applies at the location fx ^ 0. If 
the physical conditions are not met, then we generate a new 
trial point. If they are met, then we proceed to test if the 
analytic criteria A to G are satisfied for this accepted point. 

In total, we repeat this process until 10^ accepted points 
are found, requiring ~ 10^ simulations in total. Although the 
number of simulations may appear modest, we note that at 
each realization we must numerically compute the functions 
l{g,) and dl{p.)/dfx at 10 ^ locations, which takes substantial 
computational overhead (~ 30s per simulation). 

We find that 95.3% of the accepted points also satisfy 
criteria A to G, or a completeness of >95%. This indicates 
that our seven criteria are slightly overly-conservative, crop¬ 
ping ~ 5% of the physically permissible LDCs. If we use the 
unmodified versions of criteria F and G (see §2.6 and §2.7) 
and repeat the exercise, we find that 100 % of the physically 
valid points satisfy the criteria. However, as discussed in §5, 
this now yields an asymmetric loci of allowed LDCs, imped¬ 
ing efforts to find an efficient sampling algorithm. 


4.3 Validity Tests 

In an analogous approach to the previous tests, we begin 
by drawing uniform random samples in {C 2 ,C 3 ,C 4 } as be¬ 
fore, except that we know constrain the cuboid to the spe¬ 
cific bounds derived in §3.2 (including C 4 >0). We then 
test whether the seven analytic criteria are satisfied or not 
and if so consider the point to be accepted. We continue 
until 10 ^ accepted points are found, which we found re¬ 
quired 101,163,869 trials. Since we used the tightest bound¬ 
ing cuboid possible here, this reveals that the most efficient 
sampling possible without transforming the {c 2 ,C 3 ,C 4 } LDCs 
would be just under 1%. Whilst one could proceed in this 
way, applying an acceptance/rejection test at each realiza¬ 
tion, uniform sampling would reject over 99% of realizations, 
making such an approach highly inefficient. 

We next test whether each of these accepted points sat¬ 
isfies the physical conditions I, II and III. As before, this is 
done by evaluating the functions I{g) and dl{p)/dfl at 10 ^ 
evenly spaced locations. 

From these tests, we find that 100% of the 10^ accepted 
points satisfy the physical criteria. Therefore, the seven cri¬ 
teria are fully valid and drawing a point which satisfies them 
is guaranteed to always satisfy the physical conditions I, II 

and III. 

These tests therefore confirm that our criteria have a 
very high completeness and perfect validity. We therefore 
proceed with confidence that they provide a suitable set 
of constraints to evaluate the range of physically plausible 
LDCs. 


5 TRANSFORMATIVE GEOMETRY 
5.1 Overview 

When plotted in the original c = {c 2 ,C 3 ,C 4 } parameter space 
(Fig. 4), the loci of physically allowed LDCs resemble a ro¬ 
tated, tilted ellipse of thin but hnite width with a symmetric 
protrusion running along the semimajor axis. This complex 
morphology cannot be easily sampled from and if one wished 
to draw a set of LDCs from a uniform prior, there would be 



Figure 4. 3D plot of the allowed LDCs in the original parameter 
space: {c 2 ,C 3 ,C 4 }. All of the plotted points satisfy the physical 
conditions I, II and III, as well as the unmodified criteria A-G. 
The green points also satisfy the modified criteria A-G (whereas 
the red do not), chosen to yield a more symmetric loci of points 
and comprising > 95% of the volume. 


no alternative except to draw from the full cuboid and per¬ 
form an acceptance/rejection test using our seven analytic 
criteria. As evident from Fig. 4, the volume of allowed points 
is much less than the bounding cuboid volume, meaning that 
such a procedure would be highly inefficient. Numerical tests 
reveal that the efficiency of such a procedure is just under 
1 %, making any algorithms using this method very wasteful. 

We are therefore motivated to try and transform the 
geometry of the accepted loci of points into a more regular 
shape that we can sample from efficiently. We previously did 
this in the 2D case of quadratic limb darkening in Kipping 
(2013), but the transformative geometry required here is not 
only more complex but also includes an extra dimension. 


5.2 Rescaled LDCs 

We begin by noting that the envelope functions shown in 
Fig. 1-3 provide an excellent description for the bounding re¬ 
gion of allowed LDCs. Despite some small exceptions, we are 
motivated to exclusively use these simple envelopes rather 
than the full criteria since i) they provide a nearly perfect 
description of the loci ii) the envelopes are symmetric func¬ 
tions derived from quadratic forms iii) all of the envelopes 
come from criteria A and G alone. In practice, criterion F is 
also necessary to remove a duplicate set of solutions. 

We therefore proceed to only consider the region con¬ 
tained by criteria A, F and G, which we denote as the sim¬ 
plified region. This simplification means that the bounding 
cuboid is slightly modified to: 
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0 < C2 < 9, 

—8--6a/2 < C3 < —8 + 63/2, 

0 < C4 < 9. (42) 

As a first transformation, we re-scale the axes into a 
unitary cube by the use of di terms, defined as: 

do = C 2 / 9 , 

^ 6 V 2 - 8 -C 3 

12^2 ’ 

d4 = C4/9. (43) 


5.3 Righting the Allowed Region 

We next note that the loci (in d parameter space) resem¬ 
bles an elliptic thin disk with the semimajor axis pointing 
along the unit vector {1,1,1}. Further, the disk appears ro¬ 
tated about this unit vector, with respect to the axes of the 
reference frame. We decided to try to right the volume by 
performing a clockwise rotation of (;r/3) radians about the 
unit vector (Mi). We then perform a further rotation which 
relocates the unit vector {1,1.1} to {1,1,0} (M 2 ), followed 
by a third rotation relocating {1,1,0} to {1,0,0} (M 3 ). The 
total rotation matrix applied is described by: 


e = M 3 .M 2 .Mi.d, 

e = M.d, (44) 

where 


M = 



\/2 

1 


\ V 6 


3/3 

0 

+2 


75/ 


(45) 


Applying these transformations gives a new-coordinate 
system of: 


62 


63 


64 


36 — 24\/2 + 8c2 — 3\/2c3 + 8C4 
72^3 ’ 

C4-C2 

972 ’ 

24 — 183/2 + 2 \ f 7 . C 2 + 3c3 + 23/2C4 
363/3 


(46) 

(47) 

(48) 


After visually inspecting the loci in this transformed 
parameter space (as shown in Fig. 5), we note that the al¬ 
lowed region now resembles a cone. Motivated by this ob¬ 
servation, we proceed to transform this shape into a cone 
of symmetric proportions and with principal axes aligned to 
the transformed frame. 


5.4 The Conal Region 



Figure 5. Same as Fig. 4, except the coordinates have been 
transformed from {C 2 ,C 3 ,C 4 } —> {c 2 ,r 3 ,C 4 }. In this parametrization, 
we observe that the volume of green points resembles a cone. 


The extrema of £ 3 ( 62 , 64 ) occur when (64 — 62 ) is max¬ 
imized/minimized, as evident from Equation 48. This can 
be considered further by studying our earlier illustration 
in Fig. 2. This can be found by maximizing the function 
( 6 ^ — 62 ) with respect to 62 , which reveals the extrema is 
—3 < (64 — 62 ) < +3. We are therefore able to show that 
— < 63 < This can also be achieved by maximiz¬ 

ing/minimizing the 63 expression using the additional con¬ 
straints of criteria A, F and G. Repeating for the other terms 
we define a new boundary box of 


3V3-2V6 2(3 + 72 ) 

<62 < 


18 

1 

372 

372-4 

673 


<63 < 


<64 < 


373 

1 

372 

372-4 

673 


(49) 

(50) 

(51) 


We choose to re-scale the e-parameter space into a unit 
vector cube via: 


/2 


h 


/4 


3+3-23/6 
^ _18 

3/2 , +3 ’ 

3 / 3 + 2 


3^3 

7?’ 


64 + 


33 / 2-4 

63/3 


63 / 2-8 

63/3 


(52) 

(53) 

(54) 


As with the c parameter space, our first step is to re-scale 
the axis into a cuboid with unit lengths, requiring us to first 
compute the extrema in e space. 


At this point, we now choose to rotate the conic section 
by njA radians in a clockwise sense around the / 3 -axis, so 
that the cone’s apex is located at the origin and the cone 
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Figure 6 . Same as Fig. 4, except the coordinates have been 
transformed from {c 2 , 03 , 04 } —> {g 2 ,g 3 ,g 4 }- In this parametrization, 
the green points are well described by a cone of radius, R = 1/2 
and height, H = (—4+ 10\/2/3). 


5.5 Sampling from the Conal Region 

The samples shown in Fig. 6 appear consistent with points 
uniformly drawn from within the volume of a cone. We here 
describe the mathematical formalism by which one can com¬ 
pute such samples. 

Samples may be drawn from a cone by first consider¬ 
ing how to draw samples uniformly from within a circle. 
This well-known problem can be tackled by using polar co¬ 
ordinates and drawing a random polar angle 6 in the range 
0 -2;rrad and a random radius r from a triangular distribu¬ 
tion between 0 and p, where p is the full radius. We now 
note that the radius varies as a function of height, h, along 
the cone, such that p{h) =Rh/H. Finally, h is drawn from 
a quadratic power-law distribution from 0 to R (the full 
height), since the area of a circle increases as p^. Drawing 
a random uniform variate for ag, and between 0 and 
1 , the polar angles, height and radius of a point uniformly 
drawn from within the cone may be expressed as 


0 = Inae, 


h = Ha, 


1/3 

‘h > 


RIlyT^ 

~H 


(64) 

(65) 

( 66 ) 


Converting these into Cartesian elements, we have 


points along the e 3 -axis. We accomplish this using an addi¬ 
tional change of variables: 


_/2-/4 
g2 2 ’ 

83 =/3i 

/2+/4 


84 = 


V2 ’ 


(55) 

(56) 

(57) 


where we have additionally normalized §2 by a factor of 
3/2 to allow the loci to be symmetric on the g 2 -g 3 plane. 

In this frame, our cone now has an apex at zero, with a 
height of H = {— 4) and a radius of = 1 /V2, as shown 
in Fig. 6. Writing out the g terms relative to the original c, 
coefficients, we have 


g 2 = rsin0, 
g 3 =rcos0, 
g4 = h- 

Or more explicitly: 


(67) 

( 68 ) 

(69) 


g2 =f?a,y^ay^sin(2;rae), (70) 

gi=Ra]J^alJ^co3{27tae), (71) 

84 =Ha]J\ (72) 

We may also express the c/ coefficients in terms of the 
uniform random variates, a,-: 


1/3 

C2=^( 28(9-5^) 


83 = 


12^/2 \ 
^(C4-C2), 


For which the inverse relations are 

C2 = (^ + 5V2)g2-3g3 + (l + ^), 

C 3 = (^-1472)g2 + (2-^) 


84, 

15 




(63) 


(58) 

-l-3a//^ 

(59) 



9 ( 

(60) 

\ 

-l-3a//^i 

(61) 


(62) 

-l-3a//^ 


—6cos(2;rae)-f (3-f 103/2sin(2;rae)^ j, (73) 
632-b 396^ 

-l-3a//^(4-2lV2)sin(2;rae) ), (74) 


+ 3a!f^(6cos{27tae) + {3 + 10V2sm{27tae)j j. (75) 


One may now work in the a parameter space, draw¬ 
ing samples from within a unit cube (see Fig. 7) and then 
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Figure 7. Same as Fig. 4, except the coordinates have been 
transformed from { 02 , 03 , 04 } —> {a*,In this parametriza- 
tion, the green points are nearly uniformly distributed within a 
unit cube. One may therefore uniformly sample from the cube in 
a-space and then transform back to c-space to efficiently sample 
physical LDCs. 


converting into a physically plausible set of LDCs using the 
above. 

A unique set of inverse relations can be defined by use 
of an arc tangent accounting for the quadrant of the radical 
and the use of a floor function. We have written python 
and FORTRAN code, called ldc3, to perform these functions, 
which is publicly available at this http. 


5.6 Testing Samples Drawn from the Conal 
Region 



We have provided no formal proof that the loci of points 
in g-space is bound by a cone, nor do we explicitly claim 
so. We merely observe that the morphology of the loci most 
closely resembles this shape, from which it is possible to eas¬ 
ily draw uniform samples. We here provide two simple tests 
demonstrating that sampling from the conal region provides 
an excellent set of physical LDCs. 

First, we generated 10^ uniform random points from 
within the cone and tested whether they satisfied the phys¬ 
ical conditions I, II and III. Using R= 1/2 and H = (—4-f 
lOv^/3), we find that 97.3% of the conal region is physi¬ 
cal, or equivalently, points sampled from this region have a 
validity of 97.3%. Similarly, using the sample of 10^ valid 
points generated earlier in §4.3, we find a completeness of 
94.4%. Therefore, points samples from the conal region crop 
~ 5% of the allowed parameter space. 

Aside from validity and completeness, we also consider 
the distribution of LDCs generated from sampling the conal 
region. We find that uniform samples from the g-cone (or 
equivalently uniform points in the a-cube) yield {c 2 , 03 , 04 } 
LDCs closely matching the distribution which would re- 


Figure 8. Dashed, purple lines depict smoothed histograms of 
10® LDCs generated by uniformly sampling from the conal re¬ 
gion and then transforming from a —> c parameter space. Grey 
histograms depict 10® LDCs drawn from a random uniform prior 
in c parameter space which also satisfy the seven analytic crite¬ 
ria. The close agreement demonstrates the effectiveness of directly 
sampling from the conal region to draw physical LDCs. 


suit from uniform sampling in c-space with a simple accep¬ 
tance/rejection of criteria A-G. This is evident in Fig. 8 , 
where we compare the nearly identical distributions from 
these two approaches. One can also see in Fig. 7, that the 
valid samples plotted in a-space (which were initially drawn 
uniformly in c-space) provide an approximately uniform set 
of points within the unit cube. 
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6 DISCUSSION 

In this work, we have presented a set of seven analytic 
criteria which may be used to assess the physicality of 
LDCs associated with the Sing et al. (2009) three-parameter 
non-linear limb-darkening law. We imposed simple condi¬ 
tions that the flux is everywhere positive, monotonically de¬ 
creases from centre to limb and has a negative curl at the 
limb. Through numerical testing, we have shown that points 
naively sampled with a simple accept/reject algorithm ap¬ 
plied to our criteria are always physically valid. Addition¬ 
ally, over 95% of the physically allowed loci of LDCs (found 
through brute force numerical exploration) satisfy the seven 
criteria, demonstrating a very high completeness. Using an 
unmodified set of criteria which retains the asymmetries 
present in the loci of allowed LDCs, the completeness is 
100 %. 

Armed with these criteria, we have re-parametrized the 
LDCs such that the loci of allowed points morphologically 
resemble a regular geometric shape, specifically a cone. We 
have shown that uniformly sampling points from the conal 
region in the re-parametrized space yields physically plausi¬ 
ble and uniformly distributed LDCs to high accuracy. Specif¬ 
ically, we find a validity of 97.4% and a completeness of 
94.4%. 

Sampling from the conal region may be achieved by 
drawing a uniform random variate in the transformed space 
{o!/,,0!r,0!e}, for which the relational expressions to the orig¬ 
inal {c 2 ,C 3 ,C 4 } LDCs are provided in Equation 75. We also 
provide public code (ldc3) in python and fortran to per¬ 
form both the forward and inverse calculation between the 
parametrizations (this http). 

Our work provides, for the first time, a practical and 
efficient framework for fitting astronomical data affected by 
limb-darkening with a law supporting three degrees of free¬ 
dom. Until now, one had to limit oneself to efficient sampling 
of a two-parameter limb-darkening law (Kipping 2013) and 
go without the major improvement in accuracy provided by 
a three-parameter law, such as that of Sing et al. (2009). Al¬ 
ternatively, one would have had to explore and marginalize 
over unphysical combinations of LDCs (which we estimate 
would occur for at least 99.9% of naively sampled points) or 
numerically test the physicality of each realization of LDCs 
(again with an overhead of rejecting the vast majority of 
points). In any case, we argue that our solution provides 
major advantages and enables the community to practically 
fit more complex limb-darkening profiles for the first time. 

The only published grids of theoretical LDCs using the 
Sing et al. (2009) law comes from Sing (2010). With the Ke¬ 
pler bandpass, we find that 99.6% of the Sing (2010) tabu¬ 
lated points satisfy physical conditions I, II and III. Fur¬ 
ther more, 97.7% of these tabulated points reproduce an a 
transformed LDC within the unit cube. These values again 
demonstrate that the a parametrization can be practically 
used to explore the physically allowed LDCs. Since we now 
have an efficient strategy to fit LDCs, there is the potential 
to verify the predictions made from theoretical models by 
the study of high signal-to-noise transits in the future. 

In some applications, having “only” 97.3% of the LDCs 
being physically valid may be insufficient and one may wish 
to ensure 100% validity. Since the seven analytic criteria 
guarantee 100% validity, one may draw a set of LDCs us¬ 


ing the a parametrization and then test if this realization 
satisfies the seven criteria. In practice, criteria B and D are 
never violated by points sampled from the conal region and 
thus it is only necessary to test hve criteria. This approach 
enables a guaranteed physically plausible set of LDCs at 
minimal computational expense. To aid the community, our 
code ldc3 can perform this test (this http). 

Whilst the quadratic law (and other two parameter 
laws) will likely remain suitable for many studies, the 
analysis of high precision data increasingly demands a 
more sophisticated treatment of limb-darkening to avoid 
this issue becoming a bottleneck in obtainable accuracy 
(Epinoza & Jordan 2015). By freely fitting high-precision 
data with our a parametrization of the Sing et al. (2009) 
limb-darkening, one can have greater conhdence that the 
parameters of interest are marginalized exclusively over the 
physically plausible parameter space and limb-darkening is 
modelled in a manner more consistent with simulations from 
modern stellar atmosphere models. We also note that infor¬ 
mative priors on our a parametrization may be used as well, 
in cases where one has strong belief in the results of stel¬ 
lar atmosphere models and the star is well-characterized al¬ 
ready, or alternatively from previous posteriors derived from 
freely fitting the LDCs. For either informative or uninfor¬ 
mative sampling, the a parametrization offers an efficient 
and physically sound pathway to exploring parameter space 
when modelling limb-darkening under the three-parameter 
law. 
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